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Evolution of Massive Haloes in non-Gaussian Scenarios 
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ABSTRACT 

We have performed high-resolution cosmological N-body simulations of a concordance 
ACDM model to study the evolution of virialized, dark matter haloes in the presence of pri- 
mordial non-Gaussianity. Following a standard procedure, departures from Gaussianity are 
modeled through a quadratic Gaussian term in the primordial gravitational potential, charac- 
terized by a dimensionless non-linearity strength parameter /nl- We find that t he halo mass 
function and its redshift evolution closely follow the analytic predictions of (Matarrese et alJ 
|2Q0C)|) . The existence of precise analytic predictions makes the observation of rare, massive 
objects at large redshift an even more attractive test to detect primordial non-Gaussian features 
in the large scale structure of the universe. 
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1 INTRODUCTION 

Inflation is considered the dominant paradigm for understand- 
• 1^ ing the initial conditions for structure formation in the Universe. 
u As a consequence of the assumed flatness of the inflaton po- 
^ ■ tential, any intrinsic non-linear (hence non-Gaussian, NG) effect 
during standard single-field slow-roll inflation is generally small. 
Thus, adiabatic perturbations originated by the quantum fluctu- 
ations of the inflaton field during standard inflation are nearly 
Gaussian distributed. Despite the simplicity of the inflationary 
paradigm, however, the mechanism by which perturbations are 
generated is not yet fully established and various alternatives to 
the standard scenario have been considered. A key-point is that 
the primordial NG is model-dependent. While standard single- 
field models of slow-roll inflation lead to small departures from 
Gaussianity, non-standard scenarios for the generation of primor- 
dial perturbations in single-field or multi-fiel d inflation allow for 
a larger level of non-Gauss ianity (see, e. g. iBartolo et al] l2002l : 
iBernardeau & Uza n 2002; C hen et a02007h . Moreover, alternative 
scenarios for the generation of the cosmolog ical perturbations like 
the curvaton ( Lvth et aL 2003: Bartolo et al.|i 20041 the inhomoge- 
neous reheating ([kofmanll2003l: IPvali et all 2004h and the Dirac- 
Born-Infeld(DBI)-inflation ('Alis hahiha et al. 1200 4*) scenarios, are 
characterized by a potentially large NG level (see for a review 
iBartolo et al.ll2004l) . 

Because it directly probes the primordial perturbation 
field, the Cosmic Microwave Background (CMB) temperature 
anisotropy and polarization pattern has been considered the pref- 



erential way for detecting, or constraining, primordial NG signals, 
thereby shedding light on the physical mechanisms for perturbation 
generation. 

Alternatively, one can consider the Large-Scale Structure 
(LSS) of the Universe. This approach has both advantages and dis- 
advantages. Unlike the CMB, which refers to a 2D dataset, LSS 
carries information on the 3D primordial fluctuation fields. On 
the other hand, the late non-linear evolution introduces NG fea- 
tures on its own, that need to be disentangled from the primor- 
dial ones. This can be done in two different ways: observing the 
high-redshift universe, such as e.g. anisotropics in the 21cm back- 
ground (Pillepich et al. 2007, see also Cooray (2006)), or studying 
the statistics of rare events, such as massive clusters, which are sen- 
sitive to the tails of the distribution of primordial fluctuations. 

So far, the investigation of the effects of NG models in 
the LSS has been carried out primarily by analytic means 



Lucchin & Matarrese 1988; Kovama et al. 1999; Matarrese et al. 
2000; Robinson et al. 2000; Robinson & Baker 2000; Verde et al. 
200d I2OOII: IK omatsu e t all l2003l: IScoccimarr o et all 12004 



Amara & Refregier 2004; Ribe iro et al. 2007 : .Sefusatti et al.l2007 : 
Sefusatti & Komatsu 2007; Sad eh et al.l2007r) . N-body simulations 
with NG initial cond i tions were performed in the earl y nineties 
(Messina et"aD I199QI: iMoscardini et al. 1991; Weinber 7 & Cold 
ll992l)T but the considered NG models were rather simplistic, being 
based on fairly arbitra ry distributions for the primordial density 
fluctuations (see also iMathis et"aD |2004|) . Only very recently, 
iKang etal.l ( l2007h have investigated more realistic NG models. 
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The suite of new N-body experiments presented in this paper 
considers, for the first time, large cosmological volumes at very 
high-resolution for physically motivated NG models (as described 
in the following section). Here we focus on massive haloes and their 
redshift evolution. A more general and exhaustive presentation of 
the simulations and the main results will be done elsewhere (Grossi 
et al., in preparation). 

This paper is organized as follows. In Section [2 we introduce 
the NG model here considered. The characteristics of our N-body 
simulations are presented in Section [3] In Section |4l we compute 
the halo mass function and compare it with analytic predictions. 
We discuss the results and conclude in Section [5] 



2 NON-GAUSSIAN MODELS 

For a large class of models for the generation of the initial seeds for 
structure formation, including standard single-field and multi-field 
inflation, the curvaton and the inhomogeneous reheating scenarios, 
the level of primordial non-Gaussianity can be modeled through 
a quadratic term in the Bardeen's gauge-invariant potential 
namely 

$ = $l + /nl($l-(^l)) , (1) 

where $l is a Gaussian random field and the specific value of 
the dimensionless non-li nearity parameter / nl depends on the as- 
sumed scenario (see, e.g.jB artolo et al.|[2004 ). 

It is worth stressing that eq.([T]), even though commonly used, 
is not generally valid: detailed second-order calculations of the evo- 
lution of perturbations from the inflationary period to the present 
time show that the quadratic, non-Gaussian contribution to the 
gravitational potential should be represented as a conv olution with 
a kernel /NL(x,y) rather than a product (iBartolo et al. 20051). 
However, for |/nl| ^ 1 as assumed in this paper, all space- 
dependent contributions to /nl can be neglected and the non- 
linearity parameter can be effectively approximated by a constant. 
In this work we will take -1000 ^ /nl ^ +1000. We notice that 
owing to the smallness of the contribution of non-Gaussianity 
implied by these values of /nl is always within the percent level 
of the total primordial gravitational potential, and does not appre- 
ciably affect the linear matter power spectrum. 

Yet, this range is larger than that of —54 < /nl < 114, cur- 
rently allowed, at th e 95 per cent confidence level, by CMB data 
dSpergel et al.ll2007b . The rationale behind this choice is twofold. 
First of all the large scale structure provides observational con- 
straints which are a priori independent of the CMB. Second of 
all, /nl is not guaranteed to be scale independent, while the LSS 
and CMB probe different scales. Indeed some inflatio nary scenar- 
ios do predict large an d scale-dependent /nl (see, e.g. JChenll2005l : 
IShandera & Tv e 2006^. 

Departures from Gaussianity affect, among other things, the 
formation and evolution of dark matter haloes of mass M. One 
way of quantifying the effect is to look at the halo mass function 
which we can express as the product of the Gaussian mass function, 
TiG (M, z), times a NG correction factor, Fng (M, /nl)'- 



n(M, z, /nl) = nciM, z) Fng{M, z, /nl) • 



(2) 



^ on scales much smaller than the Hubble radius, Bardeen's gauge- 
invariant potential reduces to minus the usual peculiar gravitational poten- 
tial 



In the extended Press-Schect er scenario, the NG factor can be writ- 
ten as jMatarrese et al.ll200d hereafter MVJ) 



FivG(M, 2:,/nl) 
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where 6^ (zc) = Sc{zc) \/l — S3,m^c{zc)/S, cf\^ is the mass vari- 
ance at the mass scale M (linearly extrapolated to z = 0), and 
5c{zc) — Ac/D^(zc) with D+ the growing mode of linear den- 
sity fluctuations and Ac the linear extrapolation of the overdensity 
for spherical collapse. The quantity 5*3, m represents the normal- 
ized skewness of the primordial density field on scale M (linearly 

extrapolated to ^ = 0), namely 5*3, m = oc — /nl [see, e.g., 

eqs. (43-45) in MVJ, in which notation /nl s in Model B]; 
5m represents the density fluctuation smoothed on the mass scale 
M. The redshift dependence is only through Zc, the collapse red- 
shift, which, in the extended Press-Schechter approach, coincides 
with the considered epoch, i.e. Zc = z. 

A more robust statistics, that we will use throughout this pa- 
per, is the ratio Rng {M,z, /nl) of the NG cumulative mass func- 
tion A^ng(> z, /nl) to the corresponding Gaussian one. Since 
in the high-mass tail A^(> M, z) ^ n(M, z) x M, we can approx - 
imate i?iVG(M,z,/NL) ^ Fng{M,zJ^l) dVerde et al.l l200lh. 
We notice that, being defined as a ratio, Rng is almost indepen- 
dent of the explicit form assumed for ug^M^ z). We found, how- 
ever, that the mass function in our simulation w i th Gau ssian initial 
conditions is well fitted by the lSheth & TormenI (Il999l) relation. 

An alternative approximation for Rng can obtained start- 
ing from eq.(62) of MVJ: Rng ^ Png{> Sc\z,M)/Pg{> 
6c\z, M), where 



Png{>Sc\z,M) 



TT Jo 



(-^) 



X sin ( X6c + 



(4) 



that reduces to the Gaussian case Pg{>6c\z,M) when (Sl,) = 0. 
As we will see, this formula provides a better fit for large positive 
values of /nl, while the previous one should be preferred when 
/nl is large and negative. 



3 NUMERICAL SIMULATIONS 



We have used the GADGET-2 numerical code (ISpringelll2005h to 
perform two different sets of N-body simulations using collision- 
less, dark matter particles only. 

The first set consists of 7 different simulations characterized 
by different values of /nl but the same ACDM model with mass 
density parameter Qrn — 0.3, baryon density Q^h — 0.04, Hub- 
ble parameter h — ifo/(100 km Mpc~^) = 0.7, primordial 
power law index n = 1 and as = 0.9, in agreement with the 
WMAP first-year data dSpergel et al.ll2003h . In all experiments we 
have loaded a computational box of 500^(Mpc/h)^ with 800^ dark 
matter particles, each one with a mass of m = 2.033 X 10^° M© 
h~^. The force was computed using a softening length ei — 12.5 
h~^ kpc. To set the initial conditions we have used the same initial 
Gaussian gravitational potential $l characterized by a power-law 
spectrum P(A;) oc which was then inverse-Fourier transformed 
to get the NG term /nl (^l ~ (^l)) in real space, and then trans- 
formed back to /c- space, to account for the CDM matter transfer 
fun ction. This procedure h as been adapted from the original one 
of (iMoscardini et al.|[l99ll ) and guarantees that all the /nl mod- 
els have the same linear power spectrum as the Gaussian case, as 
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we checked in the different reaUzations. The gravitational poten- 
tial was then used to displace particles according to the Zel'dovich 
approximation, starting from a glass-like distribution. Besides the 
Gaussian case (/nl = 0), we have explored 6 different NG scenar- 
ios characterized by /nl = ±100 , ±500 and ± 1000. 

The second set of simulations is designed to match those 
of (iKang et aDl2007h . Therefore, we have used smaller boxes of 
300^(Mpc/h)^ with 128^ particles and a different choice of cos- 
mological parameters {Q.m — 0.233, r^A — 0.762, Hq - 
73 km Mp c~^ and erg = 0.74) in agreement with the WMAP 
three-year data ('Spergel et al."'2007'). We have used the same pro- 
cedure to generate the primordial NG, in or der to explore th e cases 
of /nl = (-58, 0, 134), investigated by dKang et al.ll2007h . 



4 THE HALO MASS FUNCTION 

To extract the dark matter haloes in the simulations we have used 
two standard algorithms: the friends -of -friends and the spherical 
overdensity methods. The haloes, identified using a linking length 
of 0.16 times the mean interparticle distance, contain at least 32 
particles. 

In Fig. [T] we show the spatial positions of the dark haloes (cir- 
cles) in the first set of simulations, superimposed on the mass den- 
sity field (colour-coded contour plots), in a redshift slice of thick- 
ness 25.6 Mpc/h cut across the computational box. The maps are 
shown at four different redshifts, indicated in the plots, to follow 
their relative evolution (from top to bottom). The panels in the cen- 
tral column refer to the standard Gaussian case. The two extreme 
NG cases of /nl = -1000 and /nl = +1000 are illustrated in the 
panels in the left and right columns, respectively. The mass density 
fields look very similar in the Gaussian and NG cases at all red- 
shifts. Deviations from Gaussianity become apparent only when 
considering, at a given epoch, the number and distribution of the 
top-massive, virialized objects. This is particularly evident in the 
2: = 2.13 snapshot, in which the first-forming cluster-sized haloes 
can only be seen in the /nl = +1000 and /nl = scenarios, in 
which cosmic structures are indeed expected to form earlier with 
respect to the /nl = —1000 case. 

The visual impression is corroborated by the plots of Fig. |2] 
in which we show the cumulative halo mass function of all models 
in the first set of simulations, considered at the same redshifts as 
Fig.[T] The number density of massive objects increases with /nl, 
as expected, and the differences between models increases with the 
redshift, confirming that the occurrence of massive objects at early 
epochs provides an important observational test for NG models. 

In Fig. [3] we show the logarithm of the ratio of the cu- 
mulative halo mass functions of NG and Gaussian models, 
Rng (M, ^, /nl), for our N-body simulations (symbols) compared 
to the theoretical predictions of MVJ (solid curves). Filled circles 
and triangles, that refer to /nl = +100 and —100, are compared 
to theoretical predictions of both eqs. ^ and ^ (solid and dotted 
lines, respectively). The agreement between model and simulation 
is striking. Apparent deviations are well within the Poisson errors 
that we do not show to avoid overcrowding. Models and simula- 
tions are still in agreement, within the errors, for |/nl| = 500 
(open circles and triangles) and |/nl| = 1000 (not shown in the 
figure) and out to 2; ^ 5, i.e. when the validity of the MVJ approx- 
imated expressions starts breaking down: we only notice that the 
analytic model tends to slightly overpredict the difference between 
NG and Gaussian cases. In particular, we have found that eq. ^ is 



a better approximation to large, positive /nl, while eq. (0 is to be 
preferred when /nl is negative. 

This result is at variance with that of Kans et al. (2007^ m 
which the disagreement with the MVJ predictions is already signifi- 
cant at the more moderate values of /nl = — 58 and /nl = +134, 
as it is clearly illustrated in the bottom panels of their Fig. 2. To 
investigate whether the mismatch is genuine or is to be ascribed to 
numerical effects we have repeated the same analysis using the sec- 
ond set of simulations. We still found an excellent agreement be- 
tween numerical experiments and theoretical expectations, which 
adds to the robustness of our results. 



5 DISCUSSION AND CONCLUSIONS 

The numerical experiments performed in this work show the valid- 
ity of the analytic model proposed by Matarrese et al. (2000) for 
the halo mass function in NG scenarios. The model fits the results 
of our numerical experiments over a large range of /nl and out to 
high redshifts despite the fact that non-Gaussianity is treated as per- 
turbation of the underlying Gaussian model. The good match with 
the numerical experiments is, however, not surprising, since devia- 
tions from Gaussianity in the primordial gravitational potential are 
within the percent level, despite the large values used for /nl- This 
result is in disagreement with that of Kang et al. (2007) who found 
a number density of haloes larger (smaller) than ours for NG mod- 
els with the same positive (negative) /nl values. The mismatch can 
be ascribed neither to resolution effects nor to the background cos- 
mology adopted in the numerical experiments and, perhaps, should 
be traced back to a different setup in the NG initial conditions. 

Our result has two important consequences. The good news is 
that we can use an analytical model to provide accurate predictions 
for the mean halo abundance in NG scenarios spanning a range 
of /nl well beyond the current WMAP constraints, without re- 
lying on time-consuming numerical experiments affected by shot 
noise and cosmic variance. The bad news is th at deviations from 
the Gaussian model are more modest than in the lKang et al.l (l2007h 
experiments, which makes more difficult to spot non-Gaussianity 
from the analysis of the LSS. This, however, might not be a prob- 
lem when considering the observational constraints that will be 
provided by next-generation surveys. Indeed, current hints of non- 
Gaussianity from large scale structures are rather weak as they are 
provide d by the excess power in the distributio n of 2dF galaxies at 
z - (Norberg et al."2002";"Baugh et al."2004'), for which it is dif- 
ficult to disentangle primordial non-Gaussianity from galaxy bias, 
and by the possible presence of protoclusters at z ~ 4, as traced 
by r adio gal axies surrounded by Ly-a emitters and Ly-break galax- 
ies jMilev et al. 2004), for which the measured velocity dispersion 
can be compared to theoretical predictions only including a realis- 
tic model for the velocity bias between galaxies and dark matter. 
However, the observational situation is going to improve dramati- 
cally, in many respects. Next-generation cluster surveys in X-ray, 
microwave and millimeter bands like the upcoming eROSITA X- 
ray clus ter survey (jPredehl et al.ll2006h . the South P ole Telescope 
surve y dRuhl et al .1120041) . the Dark Energy Survey dAbbott et al. 
| 2005h and the Atacama Cosmology Telescope survev dKosowskv 
TOOo) should allow us to detect clusters out to high redshift. The 
abundance of these objects is, however, more sensitive to the pres- 
ence of a dark energy component than to that of a primordial non- 
Gaussianity, at least for |/nl| ^ 100. Indeed, as pointed out by 
ISefusatti et al.l (l2007h . dark energy constraints from these surveys 
will not be substantially affected by primordial non-Gaussianity, 
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f-„=-1000 



z=5.16 



z=3.96 





Figure 1. Mass density distribution and halo positions in a slice cut across the simulation box. The color-coded contours indicate different density levels 
ranging from dark (deep blue) underdense regions to bright (yellow) high density peaks. The halo positions are indicated by open circles with size proportional 
to their masses. Left panels: NG model with /nl = —1000. Central panels: Gaussian model. Right panels NG model with /nl = +1000. The mass and 
halo distributions are shown at various epochs, characterized by increasing redshifts (from bottom to top), as indicated in the panels. 
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Figure 2. Cumulative mass function of the haloes in the N-body experiments at the same redshifts as in Fig.[T] The thick solid curve represents the halo mass 
function in the Gaussian model. Thin solid line: /nl = ±1000. Dotted line /nl = ±500. Dashed line /nl = ±100. The mass function of NG models with 
positive (negative) /nl lies above (below) the Gaussian one. Poisson errors are shown for clarity only for the Gaussian model. 



as long as deviations from the Gaussian model do not dramati- 
cally exceed current WMAP constraints. However, these surveys 
could provide a fundamental cross-check on any detection of non- 
Gaussianity from CMB, especially from Planck since the physical 
scales probed by clusters differ from that of Planck by only a factor 
of two, allowing us to detect a possible scale dependence of non- 
Gaussian features. More effective constraints on both dark energy 
and non-Gaussianity based on statistics of rare objects is likely to 
be provided by the imaging of the mass distribution from gravi- 
tational lensing of high-redshift 21 cm absorption/emission signal 
that, if observed with a resolution of a few arcsec, will allow us to 
detect haloes more massive than the Milky Way back to z ^ 10 
(lMetcalf& Wh ite 2006). 

Tight constraints on primordial non-Gaussianity can be set 
by measuring the high-order statistics of the LSS. In particu- 
lar, the bispectrum analysis of the galaxy distribution at mod- 
erate z in the next-g eneration redshift surveys such as HET - 
DEX toill et all l2004l) and WFM0S2 (iGlazebrook et alJ l2005h 
will efficient ly disentangle non-linear bia sing from primordial non- 
Gaussianity jSefusatti & Komatsull2007h . 



In this w o rk we have demonstrated the validity of the 
Matarrese et al. 1 (l200Ql) approach to detect non-Gaussianity from 
the statistics of rare objects. In a future paper we will extend our 
analysis by considering the Probability Distribution Function of 
density fluctuations, that is easily obtained from our numerical sim- 
ulations and that in principle can be determined by measuring the 
21 cm line, either in emission or in absorption, before or after the 
epoch of reionization (.Furlanetto..2006.) . 
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Figure 3. Logarithm of the ratio of the halo cumulative mass functions Rng as a function of the mass is shown in the different panels at the same redshifts as 
in Fig.[T] Circles and triangles refer to positive and negative values for /nl; open and filled symbols refer to /nl = ±500 and /nl = ±100, respectively. 
Theoretical predictions obtained starting from eqs. (|3) and (|4j are shown by dotted and solid lines, respectively. Poisson errors are shown for clarity only for 
the cases /nl = ±500. 
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